#########################################################################################
####################################  FIGURE A1   #######################################
###################   Distribution of Ambivalent Sexism by Source  ######################
#########################################################################################

library(readstata13)
library(dplyr)
library(plyr)
library(ggplot2)
library(weights)

rm(list=ls())

##Set working directory 
#setwd("~/Dropbox/StrengthInNumbersReplicationPackage/replicable/figure_a1")
#########################################################################################

sexism_uma <- read.dta13("figure_a1_dataset.dta")


##Create kdensity plots of key variables
mdat_sexism_uma <- summarise(sexism_uma, loc.mean=mean(amb_sexism, na.rm=TRUE))

mdat2_sexism_uma <- ddply(sexism_uma, "source2", summarise, loc.mean=mean(amb_sexism, na.rm=TRUE), loc.med=median(amb_sexism, na.rm=TRUE))

mdat3_sexism_uma <- ddply(sexism_uma, "source2", summarise, 
                          wmean=weighted.mean(amb_sexism, teamweight, na.rm=TRUE), 
                          uwmed=median(amb_sexism, na.rm=TRUE))


##Distribution of moral individualism by sample
png(filename = "figure_a1.png", width = 1000, height = 750)
sexism_bysource <- ggplot(sexism_uma, aes(x=amb_sexism, weight=teamweight, colour = as.factor(source2), fill = as.factor(source2))) +
  geom_density(alpha=0.6)+
  geom_vline(data=mdat3_sexism_uma, aes(xintercept=wmean,  colour=as.factor(source2)),
             linetype="dashed", size=1)+
  scale_fill_manual(values=c("#8da0cb", "#66c2a5"), labels=c("Study1+Study2", "2018 CCES"), name="Source") + 
  scale_color_manual(values=c("#8da0cb", "#66c2a5"), labels=c("Study1+Study2", "2018 CCES"), name="Source") +
  xlab("Sexism") +
  theme_minimal()
sexism_bysource
dev.off()

coefficients(summary(lm(amb_sexism ~ source2, data=sexism_uma, weight=sexism_uma$teamweight)))

wtd.t.test(x=sexism_uma$amb_sexism[sexism_uma$source2==1], y=sexism_uma$amb_sexism[sexism_uma$source2==2], 
           weight=sexism_uma$teamweight[sexism_uma$source2==1], 
           weighty=sexism_uma$teamweight[sexism_uma$source2==2],
           samedata=TRUE, bootse=TRUE)